Vaccine data modeling: Custom MCMC solution

In this notebook I’ll do all the modeling for study 1. I will use my custom scoring method with BiDAG to do MCMC search for a structure across the space of DAGs. I’ll do a few searches, imposing different top-down theory/constraints. In each case, we’ll inspect the MAP dag and the posterior probability of edges.

My approach will be to searching for an underlying cognitive model. For each family, the responses are treated as generated from a cognitive model, as reports of the credence in the states of affairs represented by each variable.

Unpacking that a bit more, the idea is that participants share a mental model of some set of states of affairs (SOAs), that is, how the states of affairs relate or go together. The relations between SOAs can be represented by conditional probability table. To estimate this CPT, I further assume these CPTS can themselves be captured by a canonical distribution–in this case, something I’m calling noisy-or-nor (needs a better name). Noisy-or-nor assumes that each parent has an independent generating or preventing influence on the child node. The probability of the child is the probability that it is generated (either by one of its generating parents, or a implied “leak” variable) and that none of its preventers prevent it. This is based on work by Lu & Yuille, and draws on a common set of assumptions in causal learning literature. The essential feature of this CPD is that it assumes independent influences of each parent on the child–this is necessary in order to estimate the parameters of the CPD from reports of the marginal probability of each SOA.

source("vacc_import_data.R", chdir = TRUE)

NAs introduced by coercionattributes are not identical across measure variables;
they will be droppedJoining, by = "ResponseId"
Setting row names on a tibble is deprecated.Warning message:
In scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  :
  EOF within quoted string

# detach("package:BiDAG", unload=TRUE)
library(tidyverse)
library(BiDAG)
# BiDAG depends on 'graph' package, source: http://www.bioconductor.org/packages/release/bioc/html/graph.html
# BiDAG depends on 'RBGL' package, source: https://bioconductor.org/packages/release/bioc/html/RBGL.html
library(bnlearn)
library(HydeNet)
library(gtools)
# source("cog-model-graph-tools.R")
source("../scripts/custom-structure-learning/cog-model-main.R", chdir = TRUE)
source("../scripts/gmod_tools.R")
source("../scripts/bnTheoryBlacklist.R")
# # customize BiDAG functions
# orig_DAGcorescore <- BiDAG:::DAGcorescore # save to switch back if desired
# orig_TableDAGscore.alias <- BiDAG:::TableDAGscore.alias # save to switch back if desired

Set up data

Participants responses were scaled from a 1 to 7 rating scale onto the open interval (0,1). We scaled these responses as recommended by [cite beta regression paper]. After scaling, we interpret these responses as indicating participants’ credence in each belief–that is, the marginal probability they assign to the state of affairs described by the scale.

# rescale beta
rescale_beta <- function(x, lower, upper) {
  # rescales onto the open interval (0,1)
  # rescales over theoretical bounds of measurement, specified by "upper" and "lower"
  N <- length(x)
  res <- (x - lower) / (upper - lower)
  res <- (res * (N - 1) + .5) / N
  return(as.vector(res))
}
d_bn_scaled <- d_bn %>%
 mutate_all(function(x){rescale_beta(x,-3,3)})
## set the seed to make your partition reproductible
set.seed(123)
trainInd <- sample(seq_len(nrow(d_bn_scaled)), size = floor(nrow(d_bn_scaled)*.80))
train <- d_bn_scaled[trainInd, ]
test <- d_bn_scaled[-trainInd, ]

Define search spaces

Initialize search space. Perform initial pruning of search space using PC algorithm. This step does not respect cognitive model assumptions. Generally, this is a heuristic step needed to get the process off the ground.

In addition to the basic level of streamlining needed, we also impose different “blacklists” to constrain the search according to theory. These different blacklists apply our “generative model” assumption with different degrees of specificity. This can (might?) help us probe where our assumptions are driving findings, versus the kinds of findings that are emerging directly out of the data.

# run to compute scoretables
# takes a few minutes if not computing scoretables ...
# takes a long time (10 hours ish) if scoretables must be computed ...
compute_scoretables <- FALSE # uncomment to recompute scoretables
source("../scripts/custom-structure-learning/fit-cog-models.R", chdir=TRUE)

[1] "skeleton ready"
[1] "score tables calculated, MCMC plus1 starts"
Time difference of 4.041459 secs
[1] "MCMC plus1 iteration 2"
Time difference of 8.41198 secs
[1] "MCMC plus1 iteration 3"
Time difference of 13.5025 secs
[1] "MCMC plus1 iteration 4"
Time difference of 18.32452 secs
[1] "skeleton ready"
[1] "score tables calculated, MCMC plus1 starts"
Time difference of 4.790354 secs
[1] "MCMC plus1 iteration 2"
Time difference of 9.972552 secs
[1] "MCMC plus1 iteration 3"
Time difference of 14.54315 secs
[1] "skeleton ready"
[1] "score tables calculated, MCMC plus1 starts"
Time difference of 4.438278 secs
[1] "MCMC plus1 iteration 2"
Time difference of 9.59913 secs
[1] "MCMC plus1 iteration 3"
Time difference of 15.01579 secs
[1] "MCMC plus1 iteration 4"
Time difference of 21.71949 secs
[1] "skeleton ready"
[1] "score tables calculated, MCMC plus1 starts"
Time difference of 5.071943 secs
[1] "MCMC plus1 iteration 2"
Time difference of 9.969839 secs
[1] "MCMC plus1 iteration 3"
Time difference of 14.65087 secs
[1] "skeleton ready"
[1] "score tables calculated, MCMC plus1 starts"
Time difference of 5.640819 secs
[1] "MCMC plus1 iteration 2"
Time difference of 11.30986 secs
[1] "MCMC plus1 iteration 3"
Time difference of 18.73676 secs

nCores <- parallel::detectCores(logical=TRUE)
# load scoretables computed by fit-cog-models.R
scoretables <- list(
  readRDS("scoretable-unconstrained.rds"),
  readRDS("scoretable-intent_is_dv.rds"),
  readRDS("scoretable-abstract_are_parents.rds"),
  readRDS("scoretable-abstract_are_parents_intent_dv.rds"),
  readRDS("scoretable-theoryBasedHierarchy.rds")
  )
# # DO MCMC - comment to skip and load saved results
# all_results <- list(list(), list(), list(), list(), list())
# 
# customize_BiDAG() # customize bidag!
# 
# for (i in 1:length(startspaces)) {
#   all_results[[i]] <- parallel_orderMCMC(train,
#                                        startspaces[[i]],
#                                        scoretables[[i]],
#                                        MAP = FALSE, # might error
#                                        1e6,
#                                        100,
#                                        blacklist = blacklists[[i]],
#                                        chains = 4,
#                                        cores = 4
#                                        )
# 
# }
# 
# reinit_BiDAG() # set it back!
# saveRDS(all_results, "allresults.rds")
all_results <- readRDS("allresults.rds")

cannot open compressed file 'allresults.rds', probable reason 'No such file or directory'Error in gzfile(file, "rb") : cannot open the connection

For all four

Now let’s do this for all five models …

plots <- list()
all_edge_df <- data.frame()
for (i in 1:5) {
  chains <- all_results[[i]]
  merged_chain <- merge_chains(chains)
  map_dag <- extract_mapdag(merged_chain)
  
  h <- HydeNetwork(as.formula(bnlearn_to_hyde_string(map_dag)))
  plots[[i]] <- plot(h)
  
  edge_posterior <- edges.posterior(merged_chain$incidence, pdag=FALSE)
  rownames(edge_posterior) <- nodes
  colnames(edge_posterior) <- nodes 
  edge_df <- reshape::melt(edge_posterior) %>%
    rename(x = X1, y = X2) %>%
    mutate(constraints = i)
  all_edge_df <- bind_rows(all_edge_df, edge_df)
}

NOTE: chains did not converge for #2, intent as DV …

print("unconstrained")

[1] "unconstrained"

plots[[1]] # intent is dv

print("intent is dv")

[1] "intent is dv"

plots[[2]] # intent is dv

print("hb + naturalism are parents")

[1] "hb + naturalism are parents"

plots[[3]] # hb + naturalism are parents

print("hb + naturalism are parents, intent is dv")

[1] "hb + naturalism are parents, intent is dv"

plots[[4]] # hb + naturalism are parents, intent is dv

print("full abstraction hierarchy")

[1] "full abstraction hierarchy"

plots[[5]] # full abstraction hierarchy

all_edge_df %>%
  mutate(constraints = factor(constraints, 
                              levels = c(1,2,3,4,5), 
                              labels = c("unconstrained",
                                "intent is dv", 
                                         "hb + nat are abstract", 
                                         "intent is dv\nand hb+nat are parents", 
                                         "theory based"))
         ) %>%
  ggplot(aes(x=x, y=y, fill=value)) +
  geom_tile() +
  facet_wrap(~constraints) +
  scale_fill_viridis_c(option="C") +
  theme_minimal() +
  theme(aspect.ratio=1) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust=.5)) +
  labs(x = "From", y = "To", title="Edge posterior probabilities", fill="p(edge)") 

Comparing posteriors …

The code below can be used to compare the posterior edge probabilities under different constraints. It’s still a bit tricky to get a sense of, but this gives the answer more or less …

all_edge_df %>%
  spread(constraints, value) %>%
  mutate(diff = `2`-`1`) %>% # which to compare ...
  ggplot(aes(x=x, y=y, fill=diff)) +
  geom_tile() +
  # facet_wrap(~constraints) +
  scale_fill_viridis_c(option="C") +
  theme_minimal() +
  theme(aspect.ratio=1) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust=.5)) +
  labs(x = "From", y = "To", title="Edge posterior prob. diffs", fill="prob. diff") 

Impressions

workign?

Focusing in …

Using theory-based model …

Extract results and plot a traceplot. For unconstrained search, chains did not all converge

BiDAG_traceplot(all_results[[5]], hide = 1, show_chains=c(1,2,3,4))

# BiDAG_traceplot(all_results[[3]], hide = 100, show_chains=c(3,4))

plot map dag

merged_chains <- merge_chains(all_results[[5]][1:4])
map_dag <- extract_mapdag(merged_chains)
h <- HydeNetwork(as.formula(bnlearn_to_hyde_string(map_dag)))
plot(h)

plot posterior

edge_posterior <- edges.posterior(merged_chains$incidence, pdag=FALSE)
rownames(edge_posterior) <- nodes
colnames(edge_posterior) <- nodes 
edge_df <- reshape::melt(edge_posterior) %>%
  rename(x = X1, y = X2)
edge_df %>%
  ggplot(aes(x=x, y=y, fill=value)) +
  geom_tile() +
  scale_fill_viridis_c(option="C") +
  theme_minimal() +
  theme(aspect.ratio=1) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust=.5)) +
  labs(x = "From", y = "To", title="Edge posterior probabilities", fill="p(edge)") 

Predictions

Now we do one-node held-out prediction as I have in the past.

# this now takes ~20 minutes or so
model_string <- suppressWarnings(write_jags_model(map_dag, train))
cat(model_string)

model{
mu.nat = 0.44133
phi.nat = 4.48041
alpha.nat = mu.nat*phi.nat
beta.nat = (1-(mu.nat))*phi.nat
nat ~ dbeta(alpha.nat, beta.nat) T(0.000001,0.999999)
mu.hb = (1-(1-0.19009)*(1-nat*0.83732*1)) * (1-nat*0.83732*(1-1))
phi.hb = 3.86845
alpha.hb = mu.hb*phi.hb
beta.hb = (1-(mu.hb))*phi.hb
hb ~ dbeta(alpha.hb, beta.hb) T(0.000001,0.999999)
mu.overpar = (1-(1-0.40728)*(1-nat*0.42131*1)) * (1-nat*0.42131*(1-1))
phi.overpar = 4.92067
alpha.overpar = mu.overpar*phi.overpar
beta.overpar = (1-(mu.overpar))*phi.overpar
overpar ~ dbeta(alpha.overpar, beta.overpar) T(0.000001,0.999999)
mu.infantImmWeak = (1-(1-0.36173)*(1-overpar*0.89569*1)*(1-hb*0.19579*0)) * (1-overpar*0.89569*(1-1))*(1-hb*0.19579*(1-0))
phi.infantImmWeak = 3.82072
alpha.infantImmWeak = mu.infantImmWeak*phi.infantImmWeak
beta.infantImmWeak = (1-(mu.infantImmWeak))*phi.infantImmWeak
infantImmWeak ~ dbeta(alpha.infantImmWeak, beta.infantImmWeak) T(0.000001,0.999999)
mu.parentExpert = (1-(1-0.04619)*(1-nat*0.67305*1)*(1-hb*0.22222*1)) * (1-nat*0.67305*(1-1))*(1-hb*0.22222*(1-1))
phi.parentExpert = 4.87873
alpha.parentExpert = mu.parentExpert*phi.parentExpert
beta.parentExpert = (1-(mu.parentExpert))*phi.parentExpert
parentExpert ~ dbeta(alpha.parentExpert, beta.parentExpert) T(0.000001,0.999999)
mu.medSkept = (1-(1-0.13495)*(1-parentExpert*0.46591*1)*(1-nat*0.69152*1)) * (1-parentExpert*0.46591*(1-1))*(1-nat*0.69152*(1-1))
phi.medSkept = 6.99253
alpha.medSkept = mu.medSkept*phi.medSkept
beta.medSkept = (1-(mu.medSkept))*phi.medSkept
medSkept ~ dbeta(alpha.medSkept, beta.medSkept) T(0.000001,0.999999)
mu.vaccEff = (1-(1-0.98011)*(1-parentExpert*0.41777*0)*(1-nat*0.29718*0)) * (1-parentExpert*0.41777*(1-0))*(1-nat*0.29718*(1-0))
phi.vaccEff = 3.08906
alpha.vaccEff = mu.vaccEff*phi.vaccEff
beta.vaccEff = (1-(mu.vaccEff))*phi.vaccEff
vaccEff ~ dbeta(alpha.vaccEff, beta.vaccEff) T(0.000001,0.999999)
mu.diseaseSevere = (1-(1-0.25372)*(1-vaccEff*0.84688*1)*(1-infantImmWeak*0.67695*1)) * (1-vaccEff*0.84688*(1-1))*(1-infantImmWeak*0.67695*(1-1))
phi.diseaseSevere = 4.00263
alpha.diseaseSevere = mu.diseaseSevere*phi.diseaseSevere
beta.diseaseSevere = (1-(mu.diseaseSevere))*phi.diseaseSevere
diseaseSevere ~ dbeta(alpha.diseaseSevere, beta.diseaseSevere) T(0.000001,0.999999)
mu.vaccDanger = (1-(1-0.04665)*(1-vaccEff*0.8508*0)*(1-parentExpert*0.84782*1)*(1-nat*0.68359*1)*(1-medSkept*0.92671*1)) * (1-vaccEff*0.8508*(1-0))*(1-parentExpert*0.84782*(1-1))*(1-nat*0.68359*(1-1))*(1-medSkept*0.92671*(1-1))
phi.vaccDanger = 6.34252
alpha.vaccDanger = mu.vaccDanger*phi.vaccDanger
beta.vaccDanger = (1-(mu.vaccDanger))*phi.vaccDanger
vaccDanger ~ dbeta(alpha.vaccDanger, beta.vaccDanger) T(0.000001,0.999999)
mu.diseaseRare = (1-(1-0.73685)*(1-diseaseSevere*0.71628*0)) * (1-diseaseSevere*0.71628*(1-0))
phi.diseaseRare = 4.67457
alpha.diseaseRare = mu.diseaseRare*phi.diseaseRare
beta.diseaseRare = (1-(mu.diseaseRare))*phi.diseaseRare
diseaseRare ~ dbeta(alpha.diseaseRare, beta.diseaseRare) T(0.000001,0.999999)
mu.vaccIntent = (1-(1-0.07796)*(1-vaccEff*0.89912*1)*(1-vaccDanger*0.44421*0)*(1-diseaseSevere*0.77511*1)) * (1-vaccEff*0.89912*(1-1))*(1-vaccDanger*0.44421*(1-0))*(1-diseaseSevere*0.77511*(1-1))
phi.vaccIntent = 4.86298
alpha.vaccIntent = mu.vaccIntent*phi.vaccIntent
beta.vaccIntent = (1-(mu.vaccIntent))*phi.vaccIntent
vaccIntent ~ dbeta(alpha.vaccIntent, beta.vaccIntent) T(0.000001,0.999999)
mu.vaccTox = (1-(1-0.01349)*(1-vaccEff*0.1554*0)*(1-vaccDanger*0.98783*1)*(1-nat*0.18824*1)) * (1-vaccEff*0.1554*(1-0))*(1-vaccDanger*0.98783*(1-1))*(1-nat*0.18824*(1-1))
phi.vaccTox = 6.06736
alpha.vaccTox = mu.vaccTox*phi.vaccTox
beta.vaccTox = (1-(mu.vaccTox))*phi.vaccTox
vaccTox ~ dbeta(alpha.vaccTox, beta.vaccTox) T(0.000001,0.999999)
mu.vaccStrain = (1-(1-0.04601)*(1-vaccTox*0.54226*1)*(1-vaccEff*0.27933*0)*(1-vaccDanger*0.91631*1)*(1-medSkept*0.23111*1)) * (1-vaccTox*0.54226*(1-1))*(1-vaccEff*0.27933*(1-0))*(1-vaccDanger*0.91631*(1-1))*(1-medSkept*0.23111*(1-1))
phi.vaccStrain = 8.81664
alpha.vaccStrain = mu.vaccStrain*phi.vaccStrain
beta.vaccStrain = (1-(mu.vaccStrain))*phi.vaccStrain
vaccStrain ~ dbeta(alpha.vaccStrain, beta.vaccStrain) T(0.000001,0.999999)
mu.infantImmLimCap = (1-(1-0.01575)*(1-vaccStrain*0.27951*1)*(1-overpar*0.25236*1)*(1-infantImmWeak*0.78215*1)) * (1-vaccStrain*0.27951*(1-1))*(1-overpar*0.25236*(1-1))*(1-infantImmWeak*0.78215*(1-1))
phi.infantImmLimCap = 6.48228
alpha.infantImmLimCap = mu.infantImmLimCap*phi.infantImmLimCap
beta.infantImmLimCap = (1-(mu.infantImmLimCap))*phi.infantImmLimCap
infantImmLimCap ~ dbeta(alpha.infantImmLimCap, beta.infantImmLimCap) T(0.000001,0.999999)
 }

# takes about 20-30 min
start_time <- Sys.time()
# test2 <- test
predictions <- test %>% mutate_all(function(x){NA}) %>% mutate_all(as.numeric)
# net.fit <- bn.fit(map_dag, train)
for (targVar in names(test)) {
  tempDF <- test
  tempDF[[targVar]] <- NA
  tempDF[[targVar]] <- as.numeric(tempDF[[targVar]])
  
  # imputed <- impute(net.fit, tempDF, method="bayes-lw")
  imputed <- make_predictions(model_string, test, nodes_to_predict = targVar)
  predictions[[targVar]] <- imputed[[targVar]]
}
predictions <- predictions %>% gather(varName, predValue)
test_analyze <- test %>% gather(varName, value) %>% bind_cols(predictions) %>% select(varName, value, predValue)
print(Sys.time() - start_time)

Time difference of 25.51131 mins